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Abstract 


Consider  an  acyclic  undirected  network  G  ■  (V,E)  with  node 
set  V  and  arc  set  E  whose  arcs  are  subject  to  random  failure. 

Let  s  be  a  node  in  V  and  T  a  set  of  nodes  in  V  such  that  s  T . 
This  paper  presents  a  relatively  complete  and  comprehensive  de¬ 
scription  of  a  general  class  of  Monte  Carlo  sampling  plans  for 
estimating  g  -  g(s,T),  the  probability  that  s  is  connected  to  all 
nodes  in  T.  The  paper  also  provides  procedures  for  implementing 
these  plans.  Each  plan  uses  known  lower  and  upper  bounds  [B,A] 
on  g  to  produce  an  estimator  of  g  that  has  a  smaller  variance 
(A-g)(g-B)/K  than  one  obtains  for  crude  Monte  Carlo  sampling 
(B-0,  A-1 )  on  K  independent  replications.  The  paper  describes 
worst  case  bounds  on  sample  sizes  K,  in  terms  of  B  and  A,  for 
meeting  absolute  and  relative  error  criteria.  It  also  gives  the 
worst  case  bound  on  the  amount  of  variance  reduction  that  can  be 
expected  when  compared  with  crude  Monte  Carlo  sampl  i  ng . 

Two  plans  are  studied  in  detail  for  the  case  T  -  |t}.-iAn 
example  illustrates  the  variance  reductions  achievable  with  these 
plans.  The  paper  next  shows  how  to  assess  the  credibility  that  a 
specified  error  criterion  for  g  is  met  as  the  Monte  Carlo  experi¬ 
ment  progresses  and  then  shows  how  confidence  intervals  can  be 


computed  for  g.  Lastly,  the  paper  summarizes  the  steps  needed  to 


implement  the  proposed  technique. 
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Introduction 


Consider  an  acyclic  undirected  network  G  -  (V,E)  with  node 
set  V  and  arc  set  E.  Suppose  that  nodes  are  perfect  but  that 
arcs  are  subject  to  random  failure.  Let  s  be  a  node  in  V  and  T  a 
set  of  nodes  in  V  such  that  a  id  T.  The  purpose  of  this  paper 
is  to  present  a  relatively  complete  and  comprehensive  description 
of  a  general  class  of  Monte  Carlo  sampling  plans  for  estimating 
g(s,T),  the  probability  that  s  is  connected  to  all  nodes  in  T. 

The  paper  also  provides  procedures  for  implementing  these  plans. 
In  the  single  node  case  with  T  -  {t},  we  speak  of  s-t  connected¬ 
ness;  when  T  »  V  -  {s},  we  speak  of  network  connectedness. 

Since  it  is  well  known  that  direct  computation  of  g(s,T) 
generally  requires  exponential  time,  a  variety  of  alternative 
methods  have  been  proposed  to  approximate  g(s,T).  Notable  among 
these  proposals  are  those  based  on  bounds  and  on  Monte  Carlo  meth 
ods.  As  an  example  of  bounds,  Esary  and  Proschan  (1963,  Thm. 

.  1  )  describe  a  lower  bound  that  derives  from  a  network  Gi  - 
(V,E)  in  which  all  minimal  cutsets  are  disjoint  and  describe  an 
upper  bound  that  derives  from  a  network  G2  “  (V,E)  in  which  all 
minimal  paths  are  disjoint. 

Any  attempt  to  use  bounds  raises  two  important  questions. 
First,  are  the  bounds  sufficiently  tight  to  forgo  the  need  for  a 
more  precise  computation  of  g(s,T)  and  second,  how  much  time  does 
the  computation  of  the  bounds  require?  Note  that  if  the  computa¬ 
tion  of  bounds  for  a  proposed  method  takes  exponential  time  in 
either  |v|  or  |e|  in  general,  then  the  question  of  tightness  is 


merely  an  academic  one.  This  is  so  for  the  Esary  and  Proschan 
bounds . 

Monte  Carlo  sampling  is  the  other  major  approach  to  approxi¬ 
mating  g(s,T).  Whereas  all  techniques  of  evaluation  accumulate 
numerical  error  as  the  size  of  the  network  grows,  the  Monte  Carlo 
method  additionally  introduces  sampling  error.  To  appreciate  the 
significance  of  this  error,  we  note  that,  if  each  arc  is  either 
on  or  off,  there  are  2  possible  states  of  the  network.  There¬ 
fore,  K  replications  of  a  Monte  Carlo  experiment  can  at  best 
account  for  only  a  small  number  of  these  states  when  |e|  is  even 
moderate  in  size,  and  this  limitation  introduces  sampling  error. 

Since  sampling  error  for  Independent  trials  or  replications 
decreases  as  it  is  clear  that  more  replications  means  a 

smaller  error.  Therefore,  the  cost  of  a  Monte  Carlo  depends  on  K 
which  in  turn  depends  on  the  level  of  sampling  error  one  is  pre¬ 
pared  to  accept  and  on  the  cost  per  replication.  Two  concepts  of 
sampling  error  arise  in  practice,  one  absolute  and  the  other  rel 

ative.  Let  gj^  denote  an  unbiased  estimate  of  g  »  g(3,T)  based 

on  K  replications.  Then  to  achieve  an  absolute  error  no  larger 
than  c  >  0  with  confidence  level  greater  than  1-6,  one  needs  to 
coll ect 

K  -  min[k:  pr(  |gj^  -  g  |  >  e)  i  6]  (1) 

replications.  For  a  relative  error  no  greater  than  e  on  g  one 


needs 


K  -  mi n[ k :  pr ( 


g 


>  e)  i  6] 


(2) 


replications.  For  a  relative  error  no  greater  than  e  on  the  fail¬ 
ure  probability  1  -  g  one  needs 

Is  “  sl 

K  -  min[k:  pr  (  — -  >  e)  ^  6]  (3) 

1-g 

replications. 

A  major  consideration  in  using  the  Monte  Carlo  method  is  to 
achieve  one  or  more  of  these  error  objectives  while  keeping  K  and 
the  cost  per  replication  within  reason.  Variance  reducing  technl 
ques  described  in  Van  Slyke  and  Frank  (1972),  Kumamoto,  Tanaka 
and  Inoue  (1977),  Easton  and  Wong  (1980),  Kumamoto,  Tanaka,  Inoue 
and  Henley  (1980),  Karp  and  Luby  (1983)  and  Fishman  (1983a, 

1983b)  are  all  intended  to  make  K  smaller  than  it  otherwise  would 
be  if  crude  Monte  Carlo  sampling  were  used.  Among  these  only  the 
Karp  and  Luby  proposal  addresses  the  issue  of  how  cost  grows  as  a 
function  of  6  and  e.  In  particular,  their  method  achieves  (3)  in 
0(  |e|  m)  time  where  m  denotes  the  number  of  failure  sets  of  the 
network.  Since  this  approach  requires  the  determination  and  stor¬ 
age  of  all  m  failure  sets  prior  to  performing  the  sampling  experi 
ment  and  since  m  can  be  large  in  practice,  the  approach  has  a 
clear  limitation.  A  shortcoming  of  all  other  approaches  is  that 
none  provides  an  a  priori  estimate  of  how  large  K  needs  to  be  in 
order  to  achieve  a  specified  error  bound. 


The  present  paper  introduces  the  general  concept  of  a  Monte 
Carlo  sampling  plan  that  relies  explicitly  on  the  use  of  a  priori 
bounds  0  <  B  S  g(s,T)  S  A  <  1  to  estimate  g(3,T)  unbiasedly. 
Section  1  describes  the  plan  and  shows  that  every  such  sampling 

plan  leads  to  a  smaller  variance  var  g„  and  to  smaller  coeffi- 

K 

dents  of  variation  T(gj^;  A,B)  -  (var  gj^)^'^^/g  and  Y(1-gj^;  A,B)  « 

-  1/2 

(var  gj^)  /(I  -  g)  than  obtain  for  crude  Monte  Carlo  sampling. 

Moreover,  worst  case  bounds  can  be  computed  for  all  these  quanti¬ 
ties  in  terms  of  A  and  B.  These  worst  case  bounds  enable  one  to 
compute  upper  bounds  on  K  for  criteria  (1),  (2)  and  (3)  thus  pro¬ 
viding  valuable  information  prior  to  beginning  the  sampling  exper¬ 
iment.  The  paper  next  derives  the  worst  case  lower  bound  on  the 
variance  reduction  achievable  when  compared  with  crude  Monte 
Carlo  sampling. 

Each  sampling  plan  has  three  major  time  related  cost  consid¬ 
erations:  1)  the  time  to  compute  A  and  B  prior  to  sampling,  2) 

the  time  to  sample  the  status  of  each  arc  in  E  on  each  replica¬ 
tion  and  3)  the  time  to  determine  connectedness  on  each  replica¬ 
tion.  With  regard  to  this  last  cost,  we  assume  that  a  depth- 
first  search  algorithm  as  in  Aho,  Hopcroft  and  Ullman  (197^*  p. 
176)  is  used  which  takes  0(max(  |vj  ,  [e]  ))  time. 

Section  2  shows  that  these  general  sampling  plana  include 
the  Van  Slyke  and  Frank  (1972)  bounds  as  a  special  case.  When 
failures  occur  independently  with  identical  failure  probabili¬ 
ties,  the  computation  of  A  and  B  takes  0(  |e(  )  time  and  sampling 


per  replication  also  takes  0(  |e|  )  time.  Section  3  shows  how  the 
method  of  Kumamoto,  Tanaka  and  Inoue  (1977)  also  fits  into  this 
class  of  sampling  plans.  Since  their  suggested  Implementation 
can  require  exponential  time  to  compute  A  and  B  and  more  than 
0(  |e|  )  time  per  replication  for  sampling,  we  describe  an  alterna 
tive  approach  that  takes  0(  (e|  )  time  for  computing  A  and  B  and 
for  sampling  per  replication.  For  fixed  6  and  e  the  results  for 
bounds  and  timings  lead  to  (1)  in  0(max(  |v|  ,  |e|  )(A-B)^)  time,  to 
(2)  in  0(max(  |v|  ,  |e(  )(A-B)^/AB)  time  and  to  (3)  in 
0(max(  |vl  ,  IeI  )(A-B)2/(i-a)(1-B))  time. 

Section  H  illustrates  how  the  proposed  sampling  plan  works 
for  estimating  s-t  connectedness  for  a  30  arc  network  for  the  Van 
Slyke  and  Frank  (VSF)  bounds  and  for  the  Kumamoto,  Tanaka  and 
Inoue  (KTI)  bounds.  Section  5  then  shows  how  to  assess  the  credi 
bility  that  a  specified  error  criterion  for  1-g  is  met  as  K  in¬ 
creases.  Section  6  next  shows  how  one  can  derive  confidence  in¬ 
tervals  for  g  by  using  exact  sampling  theory,  Chebyschev- 1 i ke 
bounds  and  a  normal  approximation.  Section  7  summarizes  all  the 
steps  needed  to  implement  a  Monte  Carlo  sampling  plan  using  the 
KTI  bounds. 

We  begin  with  several  useful  definitions.  For  each  ieE  let 
xj  «  1  if  arc  i  operates 

-  0  otherwise 

X  -  (xj,  ieE)  -  a  state  of  the  network 
X<  -  set  of  all  states  x  such  that 


J  -  0.1,...,  |e| 


9 


p . 

p 


■ 


I  X.  2  J 

IeE 


Pi 


probability  that  arc  i  operates 


and 

P(x)  -  probability  that  state  x  occurs. 

Observe  that  if*  failures  occur  independently  then 

P(x)  -  n  [1  -  p.  +  x.(2p.-1)]  xeX  =  X 
IeE  ^  ^  ^ 


(M 


Lastly,  we  define  the  structure  function 

i(i(x)  =  (p(x;  s,T)  -  1  if  node  s  is  connected  to  all  nodes 

i  n  T 

=  0  otherwise 

and  the  s-T  connectedness  probability 

g  =  g(3.T)  »  I  (p(x)  P(x). 
xeX 


Note  that  if  T  =  |t}  then  the  existence  of  an  operating  path  from 
s  to  t  implies  connectedness.  If  T  -  V  -  (sj,  then  the  existence 
of  a  spanning  tree  implies  connectedness. 


1 .  Sampling  Plans 

Let  r  denote  the  set  of  all  sampling  plans  j  iji  (  x  )  ,  Q  ( x  ) ;  xeX} 
where  each  {\p(x)}  is  a  binary  function  of  the  form 

\^(x)  -  (A~B)(p(x)  *  B  0<BSgSA<i  (5) 

and  {q(x)}  is  a  sampling  distribution  such  that 
i(;(x)Q(x)  -  g, 


or  equivalently, 


I  ^(x)Q(x)  -  (g-B)/(A-B). 
xeX 

If  one  draws  K  independent  samples  x^  ^  ^  ,  .  .  .  ,  x  ^  from  {Q(x)}  then 

1  (J) 

Br  -  K  I  '('(x  )  (6) 

j-1 


is  an  unbiased  estimator  of  g  with 

var  gj^  -  (A-g)(g-B)/K  <  (A-B)^/'JK  (7a) 


and  coefficients  of  variation 

1/2 


Y(gR;  A.B)  - 


(var  g^) 


g 


<  (A-B)/2(KAB) 


1/2 


(7b) 


and 


XI-Sr;  A.B) 


(var  gj^) 


1/2 


1-g 


^  (A-B)/2[K(1-A) (1-B)  ] 


1/2 


(7c) 


These  results  which  follow  from  maximization  of  (A-g)(g-B), 
(A-g) (g-B)/g2  and  ( A- g ) ( g- B ) / ( 1 - g ) ^ ,  respectively,  have  several 
notable  features.  Observe  that  they  apply  for  dependent  as  well 
as  independent  arc  failures.  Also,  observe  that  min(  |A-g|  ,  |g“B|  ) 


moreso  than  A-B  determines  the  magnitude  of  var  gj^ ,  revealing 


the  greater  benefit  of  a  single  tight  bound  as  compared  to  a 
small  interval  A-B.  For  an  absolute  error  e  as  in  (1)  one  has 
K  $  [6(6  )  (A-B)/2e  ]2  (8) 


where 


6(6)  -  mini  6:  pr[ - >  6]  ^  6). 

( var  gj^ ) 

For  a  relative  error  e  in  (2)  and  (3)  one  has,  respectively, 

K  S  [6(6) (A-B)/2e  ]^/AB  (9) 

and 

K  S  [6(6 ) (A-B)/2e  ]2/ ( 1-A) ( 1-B) .  (10) 

To  compute  the  bounds  in  (8),  (9),  and  (10)  one  needs  to 
know  j6(6),  0  <  6  <  l}.  Although  S  =  K ( g^^ B ) / ( A- B )  has  the 
binomial  distribution  with  parameters  K  and  u  »  ( g-*B)/ ( A-B)  ,  this 
is  0^  limited  value  at  this  point  since  y  is  unknown.  From 
Chebyshev's  inequality  one  has 

6(5)  -  1/6^/2^  (11) 

which  leads  to  a  conservatively  larger  upper  bound  than  is  gener¬ 
ally  required. 

Observe  that  as  K  ^  ®  the  distribution  of  (  S- Ku  )  /  [Ky  (  1  -  y  )  ]  ^ '^  ^ 
converges  to  the  standard  normal  distribution.  Since  K  needs  to  be 
large  when  c  is  small,  one  may  in  such  a  case  use 

6  2 

6(6)  -  [6:  — - tT"?  /  -  (1-6)/2]  (12) 

(2^)^'*^  -- 

in  (8),  (9)  and  (10),  but  noting  its  approximating  nature.  These 
bounds  on  K  provide  a  convenient  prior  assessment  of  worst  case 
effort  and  make  clear  the  desirability  of  striving  for  short  in¬ 
tervals  [B,A]. 

Procedure  BOUNDS  describes  the  steps  needed  to  perform  K 
independent  replications  on  a  Monte  Carlo  experiment  using  a  sam- 


pling  plan  in  r.  Observe  that  is  the  maximum  likelihood 
estimator  of  g  and  V(g„),  apart  from  the  division  by  K  -  1 
instead  of  K,  is  the  maximum  likelihood  estimator  of  var  g  . 
The  use  of  K  -  1  in  place  of  K  makes  V(gj^)  unbiased  in  small 
samples  . 


Procedure  BOUNDS 
Purpose :  To  estimate  g(s,T). 

Input;  Network  G  =  (V,E,s,T),  bounds  A  and  B,  sampling 
distribution  |q(x),  xeX}  and  sample  size  K. 

Output :  Point  estimate  g^^  of  g(s,T)  and  point  estimate  V(gj^)  of 

var  g^. 


Method; 

I.  Initialization 
Start  with  S=0. 


II.  On  each  replication  do; 

a.  Determine  state  by  sampling  x  from  |q(x)). 

b.  Check  for  s-T  connectedness;  if  connected  add  1  to  S. 

III.  Computation  of  final  estimates 

a.  =  (A-B)  S/K  +  B. 

b.  V(gj^)  =  (A-B)^(1-S/K)(S/K)/(K-1  ). 


End  of  procedure. 
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1.1  Crude  Monte  Carlo  Sampling 

The  significance  of  a  sampling  plan  1 1|)  (  x  )  ,  Q  ( x  )  }  can 
stood  best  by  first  considering  the  simplest  form  of  Mont' 
sampling.  Here  A  «  1,  B  =  0  and  Q(x)  *  P(x)  so  that 

\|;  (  X  ^  ^  ^  ) . are  the  outcomes  of  K  independent  Ben 

trials  with  parameter  g  and 

var  gj^  =  g(1-g)/K. 

If  arc  failures  are  independent  then  P(x)  has  the  form  { 
that  determining  the  status  of  each  arc  is  itself  an  indej 
Bernoulli  trial. 

Now  observe  that  g ( 1 - g ) / ( A- g  )  ( g- B )  gives  the  number  i 
cations  needed  to  achieve  the  same  variance  with  crude  Moi 
Carlo  sampling  as  one  obtains  using  the  bounds  B  and  A. 
fore  the  worst  relative  performance,  with  regard  to  varia 
a  sampling  plan  based  on  B  and  A  occurs  when  this  quantit 
minimal;  namely  at 

g*  -  1/}l+[(1-A)(1-B)/AB]^''2j 

so  that  g '  ( 1  -  g  '  )/ ( A-- g '  )  ( g' -  B )  provides  an  a  priori  lower 
the  potential  variance  reduction. 

For  the  remainder  of  the  paper  we  assume  arc  failure 
independent.  Procedure  B  describes  how  one  samples  from 
P(x)l,  using  W  -  E  as  input  in  0(  [e]  )  time. 


Procedure  B 


Purpose ;  To  sample  the  status  or  each  arc  leW. 

Input ;  W  and  IPj^t  ieW}. 

Output ;  X  ■  Ixi,  ieW}. 

Method : 

I.  For  each  ieW  do:  sample  from  U(0,1)  and  set 

*1  *  L  Pi  ^  u  j. 

End  of  procedure. 


2 .  Bounds  Based  on  Minimal  Cardinality  Paths  and  Cutsets 

Crude  Monte  Carlo  sampling  makes  no  use  of  prior  information 
regarding  the  network  under  study,  a  weakness  that  can  be  re-- 
moved,  at  least  conceptually,  with  minimal  effort.  For  the  net¬ 
work  G,  there  exists  a  set  of  arcs  M  -  M(s,T)  such  that  at  least 
arcs  must  operate  in  order  for  s-T  connectedness  to  be  possi¬ 
ble.  Also,  there  exists  a  set  of  arcs  C  -  C(s,T)  such  that  at 
least  |c|  arcs  must  fail  in  order  for  s-T  disconnectedness  to  be 
possible.  For  example,  for  T  -  |t}  M  denotes  the  s-it  path  of 
minimal  cardinality  and  C  denotes  the  s-t  cutset  of  minimal  cardi¬ 
nality.  For  T  »  V  -  (sj,  M  denotes  the  s-T  spanning  tree  of  mini¬ 
mal  cardinality  and  C  denotes  the  network  cutset  of  minimal  cardi¬ 
nal  i  ty . 

The  significance  of  M  and  C  is  that  for  every  integer  L  S  |^| 
♦(x)  -  0  for  x  t  Xl_i 


and  for  every  integer  H  S  |c| 


n13~ 


4)  (  X  )  -  1 


ror  X  t  X  1^  I  -  X  |g| 


To  incorporate  this  information  into  the  sampling  procedure,  one 
sets 


A 


-  I 

xeX, 


P(x) 


(13a) 


B  -  1  -  I  P(x) 

t^l  “H 

and  samples  x  using  the  distribution 

Q(x)  -  P(x)/(A-B)  X  e  ^|e|-h'^L-1‘ 


(13b) 


(U) 


To  benefit  from  these  bounds,  one  first  needs  to  know 
and  |c|  .  The  determination  of  M  takes  0(  |e1  )  time  and  of  C  takes 

0(  |v|  2/3  |e^  )  time  ( Papadimi  tr  i  ou  and  Steiglitz,  1  982,  Th.  9.3.  P. 

213).  Presumably  L  and  H  are  then  chosen  so  that  the  calculation 
of  A  and  B  are  computationally  feasible.  For  example,  if  p^  »  p 

ieE  then  with  ease  one  can  choose  I,  »  and  H  »  |c|  and  compute 

A  -  1  -  Fl- 1 (  N  .  P ) 

and 

B  -  ^  -  ^|E|  -H^ 


in  0(  [e]  )  time.  By  contrast,  for  unequal  one  needs  to  perform 

0(  I  (  W 

i-0 

pute  B,  tasks  that  may  be  burdensome  if  and  |!3|  are  large  and 

one  takes  L  -  ^|  and  H  -  |c(  .  When  this  is  so,  choosing  smaller 

L  and  H  reduces  the  burden,  but  also  widens  the  interval  A-B. 

In  contrast  to  the  independent  Bernoulli  sampling  for  each 
arc’s  status  in  Procedure  B,  the  present  procedure  calls  for  sam¬ 
pling  X  on  each  replication  subject  to  the  constraint 

L  S  I  X  S  |e1  -  H.  (15) 

ieE 

For  the  case  of  pj^  -  p  icE,  Procedure  BIN  can  effect  this  sam¬ 
pling,  with  n  -  |Ej  ,  0  ■  p ,  a  »  ^|  ,  b  -  [l]  and  W  -  E  as  input, 
in  0(  |e|  )  time.  In  particular,  step  Ilb  determines  the  number  of 
operating  arcs  in  0(1)  time  if  one  uses  the  outpoint  method  of 
sampling  from  a  discrete  distribution  as  described  in  Fishman  and 
Moore  (1984).  This  method  requires  the  preparation  of  two  tables 
prior  to  doing  any  sampling.  These  tables,  whose  computation 
take  0(  |e|  )  time,  are  constructed  in  such  a  way  that  subsequent 
sampling  on  each  replication  takes  constant  time,  independent  of 
|c|  ,  |e|  and  ^1  .  As  an  alternative  one  may  sample  k  in  step  lib 

using  the  alias  method  of  Walker  (1977)  as  described  in  Kronmal 
and  Peterson  (1979)  in  the  same  time. 

lEj 

Once  k  is  known,  there  are  (  y  )  ways  of  assigning  the  k  oper- 


I?  lEl 

))  steps  to  compute  A  and  0(  (  '^  ) )  steps  to  com- 
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ating  arcs  among  the  |e1  arcs.  Step  Illb  selects  the  combination 
m  that  is  to  be  used  on  this  replication  in  0(1)  time  and  step  IV 
uses  the  k-canonlcal  representation  of  the  Integer  m  (Kruskal 
1963  and  Katona  1966)  to  identify  the  arcs  that  operate  in  0(  |e1  ) 
time.  Additional  details  about  this  sampling  method  appear  in 
Fishman  (1983b). 

Van  Slyke  and  Frank  (1972)  first  proposed  using  the  bounds 
in  (13)  to  reduce  variance  in  Monte  Carlo  experiments  in  the  case 
of  Independent  arc  failures  with  p^  »  p  ieE.  Hereafter  we  call 
these  the  VSF  bounds.  An  extension  of  this  approach  is  also 
possible.  One  takes  L  >  ^1  ,  H  >  |cl  and  sets 

A  -  1  -  I  P(x)  -  I  (t.(x)  P(x)  (16a) 

N 

and 

B  -  1  I  P(x)  -  I  4,(x)p(x).  (16b) 

lE-Cl  |E-Cl  ■  ^  N  -H 

The  idea  here  is  that  one  selects  L  and  H  and  uses  the  resulting 

A,B  and  {Q(x)},  as  constrained  by  (15),  to  effect  a  greater 

variance  reduction  than  L  -  ^|  and  H  -  |c|  allow.  This  nibbling 

away  at  the  reliability  computation  has  an  increasingly  beneficial 

effect  as  min  p  -»  1  since  the  term  for  which  I  x  -  |e-C|  becomes 

i  ieE  ^ 

important  in  (16a).  It  has  a  similar  effect  when  max  p .  0  since 

i  ^ 

the  term  with  'I  x  ■  |^|  becomes  important  in  (16b).  The  feasi- 

ieE 

bility  of  using  this  method  depends  on  the  ability  to  compute 


xeX 


N  '1 


xeX, 


Procedure  BIN 


Purpose ;  To  sample  the  status  of  each  arc  i  in  W. 
Input !  n,e,  |Fi(n,0);  i  -  0,1, a,b,W  -  \ 
Output ;  X  -  {xi,  leW}. 

Method : 

I.  Initialization 

a.  For  each  ieW,  set  x^  -  0. 

II.  Determine  the  number  of  arcs  that  operate 

a.  Sample  U  from  0(0,1). 

b.  k  -  min(j:  Fj(n,0)  -  Fb-l(n,0)  >  U[Fn-,a(n 

Fb-i(n,0)];  b  S  j  i  a). 

III.  Determine  the  "number"  of  the  combination 

a.  Sample  U  from  U(0,1). 

b.  m  - 

IV.  Determine  status  of  arcs 

a.  to-kandA-m. 

b.  Until  A  -  0  do:  y  -  max[z:  (  )  ^  A]; 

U) 

y 

i  »  Wy  +  ^;  Xi  ■  1  ;  A  ■  A  ~  ID  *  u) 

c.  Until  (D  -  0  do:  i  -  w  ;  x,  -  1  :  and  id  «  id 

ID  1 


End  of  procedure. 


Since  the  con  pa- 


♦  (X)  for  each  x  e  and  x  c  X  -x  |jj 

tational  feasibility  of  A  and  B  in  (16)  are  in  question  even  for 
the  special  case  of  equal  ,  the  general  applicability 
of  this  approach  is  even  more  limited. 


3 •  Bounds  Based  on  a  Control 

The  bounds  [B,A]  in  Section  2  are  global  in  nature  and  ex¬ 
ploit  a  relatively  small  amount  of  information  about  the  network 
under  study.  Moreover,  their  application  is  severely  limited  by 
the  fact  that  only  in  the  case  of  equal  probabilities  are  the 
bounds  B  and  A  easily  calculated  and  the  arc  states  easily  sam¬ 
pled.  The  present  section  considers  more  informative  bounds  that 
rely  on  the  status  of  a  subset  of  arcs  and  shows  how  the  result¬ 
ing  local  conditionality  enables  one  again  to  achieve  a  variance 
reduction.  This  approach  is  originally  due  to  Kumamoto,  Tanaka 
and  Inoue  (1977),  but  the  presentation  here  is  considerably  dif¬ 
ferent  in  organization  and  more  comprehensive  in  character.  We 
refer  to  the  resulting  bounds  as  the  KTI  bounds. 

Let  j4ii(x)}  and  denote  binary  functions  on 

1 0 , 1  }  such  that 

((il(x)S(|,(x)S(|i2(x).  (17) 

Then  for 


i 


9 


1  8' 


g,  -  I  (x)  P(x) 
^  xeX 


one  has 

gl  ^  g  S  g2* 

Observe  that 

<ji1  (x)<pix)  -  1  if  (x)  -  1 

-  0  if  <>1  (x)  -  0 


1  .2 


and 


(>2(  X  )  41  (  X  )  -  1 


if  4i(x )  -  1 

if  <}.(x)  -  0 


so  that 


and 


4il(x)4>(x)  -  4ii(x) 

(|i2(x)4i(x)  -  4i(x). 


4  (x)^4i  (x) 

Q(x)  -  - — — -  P(x)  xeX 

®2 


and  suppose  one  samples  x  from  jQ(x)).  Then 


Let 


I  4(x)Q(x) 
xcX 


gitg. 


g2'^g/ 


( 18) 


which  suggests  that  one  take  B  -  gi  and  A  -  g2.  Then  var  i|)  ( x ) 
“  (82*g)  (g"gi  )<g(  l-'g) . 

Two  remaining  issues  concern  the  selection  of  |4»^(x)}  and 
{(l>^(x)}  and  the  method  of  sampling  from  (Q(x))  in  (18).  For 


T  -  {t},  Kumamoto,  Tanaka  and  Inoue  suggest  that  one  can  deter¬ 
mine  a  lower  bound  based  on  minimal  s-t  paths  and  an 

upper  bound  {(Ji2(x)}  based  on  minimal  s-t  cutsets.  In  particular, 


let  Pi,...,Pj  denote  I  minimal  s^t  paths  and  C'(,...,Cj,  J  minimal 
s-t  cutsets  for  the  network  G  -  (V,E).  Then  (17)  holds  for  a 
coherent  system  if  one  defines  the  bounds  as 

i|),(x)-1-  S  (l-n  x)  (19a) 

J-1  ieP. 

and 

^  (x)  -  1  -  [1  -  u  (1-x  )].  (19b) 

J»1  ieC. 

Observe  that  (19a)  and  (19b)  act  as  controls  on  the  values  that 

the  structure  function  ((|i(x)}  can  assume. 

To  implement  this  approach,  one  needs  to  determine 

P^,...,Pj  and  Ci,...,Cj,  compute  g^  and  g2  and  devise  a 

sampling  plan  for  x.  Kumamoto,  Tanaka  and  Inoue  address  only 

the  sampling  issue.  Let 

I  J 

P^)u(uC.) 
j-1  J-1 

so  that  sampling  xj  for  ieE-n  involves  a  Bernoulli  trial.  Sam-- 
pling  arcs  in  n  takes  more  care.  Kumamoto,  Tanaka  and  Inoue  sug¬ 
gest  sampling  these  arcs  sequentially.  Define  two  disjoint  arc 
sets  fli  and  {I2  such  that  n  -  ♦  ^2  de^'ine 


X,  n  p,  ) 


1 


I 

n  ( 1-  n 

J-1  leP  -u  n 
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and 

J 

4i„(n  ,fi-)  -  n  [1-  n  (1-x  )  n  d'-p,)]. 
J-1  leC  uQ  ieC  uR2 


Then 


pr  (x. 


[(J)2(0.P2"  i'"  ^  ^  ^ 


®2  ~  ®1 


and  ( 20 ) 

.^2"  ’“2"  ^  ^ 

pr(x  *1|x.,iefl.)-  -  • 

^  ♦2^“i  ’”2^  ”  4.^  (n^  ,02) 

|fi  1 1  >  0  ,  r  e  fi  2  • 

One  sees  that  the  evaluation  of  4i2(*.’)  and  ij)(i(*,*)  can  take 


J  I 

0(  I  |c  J  )  time  and  0(  I  |P,|  )  time,  respectively,  so  that  for 

j-1  ^  j-1  ^ 

J  I 

all  arcs  in  p  these  computations  take  0(  |p|  (  I  |c  J  +  I  |P  J  )) 

j-1  ^  j-1  ^ 


time . 

Without  any  further  specification,  the  computation  of  gi 

I  |Pjl  I  ICjl 

and  g2  takes  0(2  )  and  0(2  )  time,  respec¬ 


tively,  and,  in  fact,  the  time  required  to  determine  I  paths 
and  J  cutsets  remains  in  question. 

To  reduce  sampling  time  when  the  determination  o^ 

Pl,...,Pl,  Ci,...,Cj,  gi  and  g2  are  computationally  feasible 


we  suggest  the  following  more  efficient  method.  Observe  that 


Q(x)  -  Q(y,z)  -  Qi(y)Q2(z) 


where 


y 

Q,  (y) 


1  Xi  ,  i  GSl }  « 
«i»2(y (y) 


z  ■  i  Xj  ,  i  eE->n  1  , 

n  [x  (2p  -  1 )  +  i-p  ] 

ieP 


and 

Qp(z)  -  n  [x  (2p  -1 )+1-p  ], 
ieE-« 

Now  one  can  compute  the 

n  -  I  [4>2(y)  "  <^1  (y)  3  (21  ) 

y  e  Y 

entries  for  (QiCy)}  in  0(2  )  time.  provided  this  is  computation¬ 

ally  feasible,  one  can  then  use  these  entries  to  compute  tables  ''or 

IeI 

the  outpoint  method  of  Fishman  and  Moore  (198^),  also  in  0(2  ''  ) 

time,  prior  to  experimentation.  Lastly,  one  needs  to  create  a  1-1 
mapping  from  the  n  outcomes  to  their  corresponding  stored  vectors  y 
«  (xj,  icO).  Then  on  each  replication  sampling  from  this  |Qi(y)) 
occurs  in  0(  In]  )  time. 

As  the  size  of  the  network  grows  the  effectiveness  of  this  meth 
od  of  variance  reduction  can  only  be  maintained  if  n  grows  with  E. 
Therefore,  there  comes  a  point  at  which  the  tabling  method  is  no 

longer  feasible;  nor  are  the  determination  of  Pi . Pj  ,  Ci,...,Cj, 

gl  and  g2  computationally  feasible  without  additional  restric¬ 
tions.  To  solve  this  problem  we  propose  that  the  paths  Pi . Pj  be 
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chosen  edge- d  i  s  j  oi  nt  and  the  cutsets  Cii.-.iCj  be  chosen  edge-'dis- 
Joint.  This  offers  immediate  advantages. 

First,  it  is  clear  that  I  ^  (:|  ,  C  being  the  s-t  cutset  of  mini 
mal  cardinality.  In  fact,  one  can  achieve  I  «  |c|  and  determine 

Pi . Pi  in  0(1*  |e|  )  time  using  a  network  flow  algorithm  with  unit 

capacities,  as  described  in  Wagner  (  1  975,  p.  95^).  Second,  JS  |^| 
where  M  denotes  the  s-^t  path  of  minimal  cardinality.  Here  one  can 
determine  the  J  »  ^|  disjoint  cutsets  in  0(  |e|  )  time  by  beginning 

at  node  s  and  appropriately  labeling  arcs.  Third,  edge-dis jointness 

I  J 

implies  that  g  and  g  are  computable  in  0(  I  |P  |  )  and  0(  I  [C  .|  ) 

j-l  ^  j=1  ^ 


times  respectively. 

Fourth,  one  can  sample  all  arcs  in  E  in  0(  Je)  )  time. 
Procedure  Q  describes  such  an  algorithm.  In  particular,  it  uses 
the  quantities: 


j  *  1  ,  I 


ieP 


UQ  -  0 

ui®  n  (^“Pj)  j»i,...,j 

I 

A  -  1  -  n  ( 1- X . ) 


u)  =  n  ( 1  -*(»). ) 

j“i  ^ 

Ji  -  j  if  lePj  j  =  1 . 1 

I 

=  0  if  I  P. 

k-1 

and 

kj«j  ifieCj 


J 

=  0  if  I 

k»1 


Here  jj^l  and  jk^)  are  pointers  identifying  the  path  and  cut¬ 
set,  respectively,  to  which  arc  i  belongs. 


4 .  Bounds  on  Bounds 

I  one  uses  a  table  or  Procedure  Q  then  an  (e,6)  absolute 
criterion  has  0(max(  jv]  ,  |e|  )  [  6  (  6  )  (  A- B  )  /  2e  ]  ^  )  worst  case  time  from 
(7a)  and  (e,6)  relative  criteria  have  0(max(  |V  |,  1e|  )  [  B  (  6  )  (  A-B  )  /  2e  3^/ AB  ) 
and  0(max(  |vj  ,  jEj  )  [  B  (  6  )  (  A  -  B  )  /  2e  ]2/ (  1  -  A )  (1  -  B  )  )  worst  case  times  from 
(7b)  and  7c)  respectively.  While  (19)  leads  to  a  variance  reduction, 
it  is  important  to  understand  exactly  how  these  worst  case  times  behave 


as  a  function  of  Ip^,  ie£)}  and  as  a  function  of  the  size  of  the  network 
G.  For  convenience  of  exposition  we  assume  |P^j  S  |  Pjj  and 


Procedure  Q 


Purpose:  To  sample  the  status  of  each  arc  i  In  E. 

Input :  E ,{ Pi  ,  IeE},  I,  J,  (Aj>* 

jwj;  j=0,1,...|j}i  A>  <*)• 

{xi  ,  ieE}  . 


Output 

Method 


I  . 
II 


III 


Initialization:  A  •  u>  and  p  =  A. 

Sample  arcs  in  E-n 

If  Pi  =  P  ^or  all  ieE-fi,  then  use  Procedure  BIN  to 
determine  status  of  these  arcs. 

Otherwise:  For  each  arc  i  in  E-Q,  sample  U  from 

U(0,1)  and  set  Xi  =  LPi*^J  • 

Sample  arcs  in  P 

For  each  arc  i  in  n  do: 

a.  Compute  the  probability  that  arc  i  fails 


1 


^  ‘"k.  ''l  ^ 


(  A-A  .  )/  (  1- A  .  ) 

J  i  J  i 


and  q  «=  (u)-A)p  /(&“p)- 
i 


Determine  the  state  of  arc  i 

Sample  U  from  U(0,1)  and  set  Xi  -  U-q  +  l  J  . 
Update  parameters  for  the  next  arc 


j/  [  1-  /  ( l-Pj  )  ] ; 


(1-Xj)w^^  /  (1-Pi) 


1- (1-A)[1-x  A  .  /p.  ];  A  -  A  X  /p  ; 

1  J  i  1  J  i  J  i  ^ 


(1)  and  p  »  A. 


End  of  procedure  . 
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|:j1  .  note  that  \p^\  >  HI  .  I  ^  HI  •  H  ll  ^  HI 
HI  .  and  begin  with  IfCgK;  A,B)  in  (7b). 


M . 1  Boundedness  with  Respect  to  jp^i  iep} 
Observe  that  for  disjoint  paths 
I 

B  =  g  =  i-  n(i-n  Pj) 
j=1  ieP. 


and  let 

r  «=  max  n  p . 
l^j^I  iePj  ^ 


and 

N  =  number  of  paths  k  for  which  n  p .  «=  r  k  = 

ieP,  ^ 
k 


Then  one 
B  *= 


can  write  (23)  in  the  more  concise  form 
Nr  +  0  (r  ) 


where  o(w)  denotes  a  function  h(w)  such  that 
Now  observe  that  for  disjoint  cutsets 


1 im  h ( w  ) 
w-»0  w 


J 

A  -  g  -  n  [  1-  n  ( 1-p.  )  ] 
j.1  ieCj 


and  let 
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q  -  max  p .  j -1  ,  .  •  .  ,  J 

J  icCj 

J 

q  -  n  q. 
j-1  ^ 

and 

M  =  number  of  p  -  q  icC  j  - 

*J  ^  *J  sj 

Then  A  has  the  more  concise  representation 
J 

A=qnM.+o(q).  (26 

j-1  ^ 

Note  that  (25)  and  (26)  hold  regardless  of  whether  or  not 

Pi . Pj  are  disjoint  and  Ci,...,Cj  are  disjoint. 

With  regard  to  (7b)  one  has  for  a  given  network  G 


(A-B)^  q  o(q)/q  -  o(r)/q]^ 

AB  *  r  J 

[  n  M,  +  o(q)/q][N  +  o(r)/r] 

J-1  ^ 


{  27 


which  is  finite  for  all  0  S  Pj  S  1  iefi  if  and  only  if 

q  -  r .  (28 

This  condition  is  met  for  Pi  -  { q j ,  j  .  For  the 

cial  case  p^  =  p  icn  only  |Pi|  «  J  is  needed. 

We  now  turn  to  T(1-gK:  A,B)  in  (7c)  for  which  one  has 
the  alternative  representations 
_  I  _ 

B-1-qnM.  ♦o(q) 

J-1  ^ 


where 


-27 


I  _ 

q  =  n  q 
j*i  ^ 


q  -  max  ( 1-p^  )  j  -  1  ,  I 

J  UPj 


Mj  »  number  of  p^,  -  1  -  qj  on  path  j  -  1 . 1 


and 


A-l-'Nr  +  oCr) 


where 


r-  max  n  (1-p.) 
1Sj<J  ieCj  ^ 


an  d 


N 


number  of  cutsets  k  for  which 


n 

ieC 


(I-P^) 

k 


r  k  «  1  ,  .  .  .  , 


Then  for  a  given  network  G 


„  [  n  M  .  -  iJ  r/q  ♦  o(r)/q  ^  o(q)/q]^ 

(A-B)^  J-1  ^ _ 

(1-A)  M-B)  -  I_ 

[  n  M  +  o(q)/q]CN  +  o(r)/r] 
j-1  ^ 


(29) 


which  is  finite  for  all  pj^  ien  i^  and  only  if 

q  -  r.  (30) 

This  condition  is  met  for  -  h~qj;  J“1 . l}.  For  the 

case  Pi  »  p  ien  only  i|  -  I  is  needed. 

These  results  carry  considerable  practical  importance  for 
they  reveal  conditions  under  which  the  time  to  achieve  an  esti¬ 
mate  of  g  or  1-g  with  specified  relative  error  is  finite  for  a 
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given  network  G  regardless  of  the  arc  probabilities  {pj,  icE}. 
Moreover,  the  forms  (27)  and  (29)  provide  valuable  guidance  when 
choosing  paths  and  cutsets.  For  example,  choosing  paths  and  cut- 

J 

sets  to  mimimize  |  n  is  a  desirable  objective  with  regard  to 

J-1  ^ 


J 

(29),  as  is  choosing  them  to  mimimize  |  n  M  -N  |  with  regard  to  (27). 

j-1  ^ 

.  2  Boundedness  with  Respect  To  Network  Size 

When  one  turns  to  the  effect  of  network  size  on  the  complex¬ 
ity  of  a  Monte  Carlo  experiment  with  an  (e,6)  accuracy  criterion, 
one  quickly  realizes  that  the  way  in  which  the  network  grows  is 
of  crucial  importance.  In  particular,  the  effect  on  g(s,t)  is  a 
central  consideration.  Since  this  topic  deserves  considerably 
more  space  than  we  can  afford  here,  we  limit  our  comments  to  a 
special  but  interesting  case. 

By  network  growth  we  mean  an  increase  in  the  number  of  arcs. 

Recall  that  M  and  C  denote  the  s-t  path  and  s-t  cutset  of  E  o^ 
mimimal  cardinalities.  If  the  sizes  of  M  and  C  remain  constant 
as  E  grows,  then  it  is  always  possible  to  choose  a  set  of  dis¬ 
joint  paths  Pi,...,Pj  I  i  1  and  a  set  of  disjoint  cutsets 
Ci,...,Cj  j  i  1  such  that  the  resulting  bounds  B  and  A  are 

I  J 

functions  only  of  the  arc  probabilities  (pi,  ie  I  P.)  and  |p.,  ie  Y  C.), 

J-1  ^  ^  J-1  J 
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respectively,  so  that  the  analysis  in  the  previous  section  applies. 
Then  a  Monte  Carlo  experiment  using  Procedure  Q  and  a  depth-first 
search  algorithm  to  determine  connectedness  meets  the  absolute  error 
criterion  (7a)  and  the  relative  error  criteria  (7b)  and  (7c)  in 
0(max(  |vl  ,  |e|  )  )  time . 

A  question  remains  as  to  how  interesting  is  the  set  of  networks 
that  satisfy  these  restrictions  on  M  and  C.  If  the  vertex  degrees  o'" 

G  are  bounded  from  above  then  the  constraint  on  C  is  not  unreasonable. 
However,  the  restriction  on  M  may  be  more  difficult  to  justify. 
Clearly,  more  remains  to  be  said  about  bounds  on  bounds  as  E  grows. 


5 .  Example 

An  analysis  of  the  network  in  Fig.  1  illustrates  the  benefits 
and  costs  of  the  VSF  and  KTI  bounds  for  s-t  connectedness.  The  net¬ 
work  has  30  arcs,  limiting  the  feasibility  of  directly  calculating  g. 

Insert  Fig.  1  about  here. 


The  example  assumes  independent  failures  with  p^  »  p 
ieE,  p  -  ,5,  .9  and  .95,  (s}  -  1  and  {t}  -  20.  For  the  VSF 
bounds  we  have  L  -  -  5  and  H  ■  jcj  -  3.  For  the  KTI  bounds 

we  consider  two  cases.  Case  1  has 

n  -  1 1  ,2,3, 't,  9, 1  1  ,  1  8,1  9.27,28,29, 30} 

Pi  -  {l  .'*.1  1  ,19,28}  P2  -  {  3,  9,1  8,27,28  } 

Cl  -  {1,2.3}  C2  -  {28.29.30} 


and  case  2  has 


Er-n 


Pi 

P3 

C2 

Cm 


(10.13.16.20.23.26} 

(3.9.18.27.28}  P2 
(2. 7. 15. 2M. 30}  Cl 
(28.29.30}  C3 
(^.5. 6. 7. 8. 9}  C5 


(1  .5.12.21 .29} 

{1.2.3} 

(11  .I2.IM.15.17.I8} 
(19.21 .22.2M.25.27} . 


Observe  that  for  case  1  |n|  -  12  enabling  us  to  use  the  tabl 

ing  method  for  (Qi(y)}.  Although  case  2  with  |o|  -  2M  limits  the 

possibility  of  tabling.  .  P2  and  P3  are  disjoint  and  C^.  C2. 

C3.  C  and  C5  are  disjoint,  enabling  one  to  use  Procedure  Q, 

Table  1  lists  the  bounds  l-sA  and  I^B  together  with  the 


Insert  Table  1  about  here. 


bounds  on  sample  size  that  (10)  induces  for  6  -  .05  and  e  ■  .05. 
Observe  that  these  worst  case  VSF  bounds  are  of  little  value  com¬ 
pared  to  the  KTI  bounds.  Also,  note  that  the  advantage  of  case  2 
relative  to  case  1,  with  regard  to  a  bound  on  K.  increases  with  p 
One  measure  of  the  effect  of  a  variance  reducing  technique 
is  the  variance  ratio  g(1-g)/K  var  gj^ .  A  ratio  greater  than 
unity  indicates  that  the  technique  has  the  desired  effect.  A  time 
ratio  T1/T2  is  also  used  where  T^  denotes  the  mean  time  to 
collect  an  observation  using  a  crude  Monte  Carlo  sampling  plan 
and  T2  denotes  the  mean  time  to  collect  an  observation  using  the 
proposed  method.  Then  Cg(1-g)/K  var  g  ]  x  Ti/T2  gives  the 


relative  time  required  to  achieve  a  given  variance  with  crude 
Monte  Carlo  sampling  as  compared  to  the  time  required  to  achieve 
this  same  variance  with  the  proposed  method. 


Table  2  shows  these  ratios  for  the  VSF  and  KTI  bounds  and 


p  »  .5,  .9,  .95  for  K  ■  2^6  ■  2621  **'i  independent  replications 
using  Procedure  BIN  for  the  VSF  case,  a  table  of  lQi(y)}  together 
with  Procedure  BIN  for  KTI  case  1  and  Procedure  Q  for  KTI  case  2. 

Insert  Table  2  about  here. 


The  results  indicate  that  both  the  VSF  and  KTI  bounds  have  advan¬ 
tages,  but  the  KTI  bounds  clearly  dominate.  Moreover,  observe 
that  Procedure  Q  consumes  considerably  more  time  than  tabling 
does.  Also  note  the  a  priori  minimal  variance  ratio 
g' (  1-g*  )/ (  A-g'  )  (g*  iB)  in  Table  2. 


6 .  Assessing  Credibility 

Since  the  worst  case  bounds  of  Section  2  often  lead  to  con¬ 
siderably  larger  sample  sizes  K  than  are  required  in  specific 
analyses,  one  is  interested  in  ways  of  assessing  the  extent  to 
which  a  specific  criterion  is  being  met  as  the  sampling  experi¬ 
ment  progresses.  Although  a  statistical  literature  does  exist 
for  sequential  sampling  to  achieve  a  fixed  width  confidence  inter 
val  ,  (e.g.  Chow  and  Robbins  1  965)  the  fact  that  it  relies  on  asym 
ptotic  behavior  (e  -►  0)  limits  its  appeal.  As  an  alternative  we 
describe  how  one  can  assess  the  credibility  of  an  absolute  or 
relative  error  criterion  as  the  experiment  proceeds. 

Recall  that  g  -  (A-'B)y  ♦  B.  Therefore,  criteria  (1),(2)  and 
(3)  imply  y  e[y  ,p*]  where  for  (1) 


for  (2) 


.  .ax[0.^(|  -  f^)] 


and  for  (3) 


w* 

« 

u 


max  I  0 
mi  n  { 1 


1 

*  1  *c 


c (1-B) 
‘  A-B 


]} 


e  (1-B)  ,  1 
A-B  ■'  ' 


Since  S  has  the  binomial  distribution  one  has 

Pi"  1  C  P,  .  P*  3  }  -  1  -  Fs(K,)i^)  +  Fs(K,u*),  (31) 

F  being  the  binomial  distribution  function. 

Table  3  shows  this  probability  for  a  relative  error  criteri' 
on  e  -  .05  for  1-‘g  for  case  2  of  the  KTI  bounds.  By  choosing 
powers  of  2  for  K  we  insure  that  at  each  successive  evaluation 
half  the  data  provides  new  information.  Observe  that  one  may  be 


Insert  Table  3  about  here. 


confident  in  terminating  the  sampling  experiment  at  K  «=  102^4  for 
p  -  .5  and  at  K  -  32768  for  p  -  .9  and  .95.  IMSL  (1982)  and  SAS 


(1982)  provide  routines  for  evaluating  the  binomial  distribution 
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7 •  Confidence  Intervals 

While  the  credibility  analysis  of  Section  6  provides  consid¬ 
erable  guidance,  one  would  also  like  to  evaluate  the  precision  of 
gj(  more  explicitly  at  the  end  o^  an  experiment.  Confidence 
intervals  provide  one  way  of  making  this  assessment.  As  the  next 
two  subsections  show,  different  methods  exist  for  constructing 
these  intervals.  Section  7.1  describes  an  exact  method  that  re¬ 
quires  some  care  with  regard  to  computation.  Section  7.2  de¬ 
scribes  a  relatively  simple  method  whose  results  rely  on  either 
general  bounds  or  asymptotic  limits. 


7.1  Exact  Confidence  Intervals 

Suppose  that  exactly  S  successes  occur  in  K  independent  trials. 
Then  there  exist  2-tuples  (y^,  , )  and  {^2’  “2^  with  0  <  ^  ^2  ^ 

and  0  <  <  02  <  1  such  that 

1  -  F5;_^(Kiy^)  “ 

and 

Fg(K,y2)  -  1  -  02* 

Let 

0(a^  )  -  (y^  ;  (K,y^  )  -  1  -  ) 

0(02^  *■  (^2*  Fg(K,y2)  ■  1  -  O2). 


and 


(32) 


-3i»- 

Then  [0{a^),  Q(a^)]  is  a  confidence  interval  for  y  at  confidence 
level  a  -  02  ~  Oi •  Since  there  are  many  <  02  that  satisfy  a  = 
02  ~  .  one  way  to  proceed  is  to  choose  o^  and  O2  to  minimize 

0(02)  -  0(a^),  subject  to  0  <  o^  <  O2  <  1  and  0=02“  ,  thus 

giving  the  shortest  interval  at  level  o.  An  alternative  way  to 
proceed  is  to  choose  o^  -  (1-o)/2  and  cx^  =  (1+o)/2  so  that  each 
tail  has  equal  probability. 

To  facilitate  the  determination  of  0(a^)  and  0(02)  in  (32) 
observe  that  (Abramowitz  and  Stegun  1964,  p.  945) 

0  . 

1  -  F  (n,6)  -  /  (”)w^  '(1-w)"  "^dw.  (33) 

J  '  0 

Therefore  0(o^)  calls  for  an  evaluation  of  the  inverse  Beta 
distribution  with  parameters  S  and  K  ^  S  +  1  and  0(02)*  ^or  an 
evaluation  of  the  inverse  Beta  distribution  with  parameters  S  +  1 
and  K  r  s.  Procedure  INTERVAL  describes  how  to  compute  the  short¬ 
est  confidence  interval  at  level  a.  Although  standard  computing 
packages  exist  to  perform  the  inversions  for  0(a^)  i  =  1,2  in  step 
lie,  experience  with  such  packages  in  IMSL(1982)  and  SAS(1982), 
indicates  that  a  large  K  severely  limits  their  abilities  to  pro¬ 
duce  numerical  results  when  S<<K.  It  is  precisely  in  such  cases 
that  one  may  be  willing  to  settle  for  bounding  or  approximating 


res  ul ts  . 
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Procedure  INTERVAL 

Purpose :  To  compute  a  confidence  interval  for  g  at  level  a. 

Input :  Sample  size  K,  number  of  successes  S,  bounds  A  and  B, 

confidence  level  a  and  grid  size  N, 

Output :  o  confidence  interval  Cg,j.g*3  of  shortest 

length. 


Method : 

I.  Initialization 

a.  L-0andU”1. 

II.  Search  for  shortest  interval 


For 

each  i  » 

1  ,  . .  .  ,N-1 

do : 

a . 

-  (1-a 

) ( i/N  )  . 

b . 

“  a  * 

c  . 

For  J  »  1 

,2  solve 

(31  ) 

for  0  (  Uj  ) . 

d  . 

If  U  -  L 

>  9(a2)  - 

9  ( 

),  then  set  L  =  0(a^)  and 

U  «  0  (  ) 

• 

III.  Compute  confidence  interval 

a.  g^  =  (A-B)L  +  B. 

b.  g*  =  (A-B)U  B. 


End  of  procedure. 


•  • 

ii«V. 


■  . 


7.2  Intervals  Based  on  Probability  Inequalities  and  Limits 


As  Section  1  indicates,  inequalities  and  approximations 
exist  that  can  be  used  to  compute  worst  case  upper  bounds  for  K 
for  specified  6  and  e.  These  inequalities  also  are  of  value  when 
computing  confidence  intervals.  Consider  the  probability  state¬ 
ment 


pr 


B^-el 


[  (  A-g)  (g-B)/K] 


1  12 


^8(1-a))^a  0<a<1 


(  3 


By  rearranging  terms  in  the  argument  of  pr ( • )  in  (2^).  one  can 
assert  that  with  probability  of  at  least  a  the  interval 


?  2  1/2 
(A-B)S->KB-f(A>B)B/2±B(A-B)[B/^<-^S(K-S)/K] 


(3 


covers  g,  where  B  =  6(1“a).  As  in  the  case  of  worst  case  bounds 

one  determines  B(1~a)  for  a  specified  confidence  level  a  from 
(11)  for  Chebyshev's  inequality  and  from  (12)  for  the  normal 
approximation.  Table  ^  presents  the  bounding  and  approximating 
.95  confidence  intervals  for  1-g  for  the  network  in  Fig.  1. 


8 .  Essential  Steps  for  Implementation 


The  results  of  this  paper  show  that  an  effective  Monte  Carlo 
procedure,  based  on  bounds  B  and  A,  exists  for  estimating  g(s,T) 
in  worst  case  time  0(max(  (v|  ,  |e|  )  •  D )  time  where  D  -  (A-B)^  for 

a  given  (e,6)  absolute  error  criterion  and  D  » 

( A- B )  2/ mi  n  [  AB  ,  (  1  -  A  )  (  1 -'B)  ]  for  a  given  {e,6)  relative  error  criter 
ion.  Moreover,  the  number  of  replications  required  to  achieve 
these  worst  case  bounds  can  be  computed  prior  to  experimentation 
and  used  as  a  guide.  Of  the  two  methods  of  deriving  bounds  de¬ 
scribed  here,  the  example  of  Section  ^  shows  that  the  KTI  bounds 
are  most  beneficial,  especially  for  p  close  to  unity.  Therefore, 
we  recommend  that  the  KTI  bounds  be  used  in  practice  together 
with  the  following  steps: 

1.  Determine  a  set  of  edge-disjoint  s-t  paths  Pi,...,Pi. 

2.  Determine  a  set  of  edge- di s j oi nt  s-t  cutsets  Ci,...,Cj. 

I 

3.  Compute  B  from  Pi,...,Pi  in  0(  |  P,  |)  time  and  A  from 

j»1  ^ 

J 

C  ,  .  .  .  ,C  ^  in  0(  1  J]  C.|  )  time. 

1  J  j  =  1  J 

Choose  a  (e,6)  absolute  or  relative  error  criterion 
((1),(2)  or  (3))  and  assign  values  to  e  and  6.  Then 
determine  the  bound  K*  on  K  accordingly  from  (8),  (9)  or 

(10).  Use  this  value  for  guidance. 

I  J 

5.  If  |Dl  -  |(  I  P,)u(  I  C ,  )|  is  small  enough,  then  compute  a 
J-1  ^  J-1  ^ 


table  of  the  probability  mass  function  iQiCy)}  in 


0(2  •'‘I  )  time. 


6.  For  the  sampling  experiment  use  Procedure  BOUNDS, 
a  .  Table  approach : 

Sample  arcs  in  E-D  using  Procedure  B  if  the  p^  are 
distinct  and  using  procedure  BIN  i^  they  are  equal. 
Sample  arcs  in  0  ^rom  a  table  of  |Qi(y)}. 

b.  Sequential  approach:  Use  Procedure  Q. 

Steps  a  and  b  each  take  0(  |Ej  )  time. 

c.  Check  for  connectedness  using  a  de  pt  h-i i  r  s  t  search 

in  0(max(  |v|  ,  |e|  ))  time. 

7.  If  the  sampling  experiment  is  to  be  performed  in  blocks 
of  Ki,  Ki  +  K2i  +  K2  +  K3,  replications,  etc.,  then 
a^’ter  each  block  compute  the  credibility  probability 
(31).  Suggested  increments  are  Kj^  =  Ki2^”^  for  i  = 

1,2,... 

8.  After  completion  the  experiment  compute  a  confidence 
interval  for  g,  as  in  Section  6. 
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Bounds  for 

Table  1 

Relative  Accurao 
6  •  .05  and  e 

y  Criterion 

-  .05 

(10) 

1-B 

1-A 

Bounds 

Chebyshev 

on  K 

Normal 

.1000x10^ 

.9395 

.9092 

.297^x10"’^ 

.2697 

67x10^ 

i»515 

3335 

13x10^ 

867 

6^0 

.5886 

.206i* 

.6867x10“^ 

.1828x10"21 

.1999x10^2 

.2002x10"2 

6^1x1023 

202522 

6il659 

12x1023 

38898 

12i4‘'9 

.1878 

.8269x10‘’1 

.1158x10"'’ 

.3352x10"29 

.3031x10"3 

.2500x10^3 

11  XI 03"’ 
5^)1635 
88683 

21x1030 

101671 

17033 

VSF 
KTI  1 
KTI  2 

.7104 

.7094 

.7105 

1 . 

1 . 
2. 

00 

89 

35 

1.00 

1.88 

2.31 

.94 

1  .04 
.24 

.94 

1.97 

.56 

P=.90 

VSF 

.2933x10''2 

1 . 

70 

1.70 

.91 

1.55 

KTI  1 

.2862x1 0“2 

16. 

67 

5.83 

1 .00 

16.67 

KTI  2 

.2867x1  0‘"2 

50. 

22 

20.93 

.26 

13.06 

P=.95 

.2938x10"3 

VSF 

5. 

33 

5.32 

.82 

4.37 

KTI  1 

.2887x10"3 

95. 

04 

13.63 

.97 

92.19 

KTI  2 

.2940x10“3 

592. 

23 

118.46 

.23 

135.49 

t  g  is  estimated  by 

estimated  by  V(g„). 

8k 

the  quantity  g(1 

-g)  and 

var  g|^  is 

Table 


.95  Confidence 

Intervale  for  1-g 

. . 

Chebyshev 

Normal 

— 

KTI  1 

p«.50 

“ 

lower 

.7065 

.7087 

upper 

.7123 

.7107 

■ 

width 

.5775x10^2 

.2531x10'"2 
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